home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Precision Software Appli…tions Silver Collection 4
/
Precision Software Applications Silver Collection Volume 4 (1993).iso
/
stats
/
chadyn.exe
/
YU.C
< prev
next >
Wrap
Text File
|
1988-12-06
|
19KB
|
674 lines
/******************* (C) 1986,7,8 by JAMES A. YORKE **************************/
/********************************** YU.C ***********************************/
#include "yinclud.h"
static double diagonal; /* for unstable manifold routines in this file
*/
/* int goodIterate; the number of iterates before leaving the box
(set in terms of diameters[] ) */
static double ddx,
ddy;
static double PreviousAdd,
oneOverX,
oneOverY;
static int multiplier;
static double jump,
oldX,
oldY,
dif;
static int itermax;
LineIterator(iterMax, iterMultiplier)/* The images of the line from
x1_line,y1_line to x2_line,y2_line are
drawn, where iterMultiplier is a number >=1
such that the line is iterated in multiples
of iterMultiplier; in the computation of
unstable manifolds of periodic orbits, it
should be 2*period where period is the
period of the orbit; the factor of two is
because the orbit might have a negative
expanding eigenvalue; iterMax is the maximum
number of times "iter" is increased; the
actual maximum number of iterates is
iterMax*iterMultiplier; routines in this
file: set_line(),LineForB(), and LineForU()
set the points(line_xi,line_yi) for i = 1
and 2. */
int iterMax,
iterMultiplier;
{
int isInBox(), level0;
int checkFreq();
double minAdd = 1;
double getAdd();
double fracY[EQUATIONS];
double oldLength = 0.,
entropy,
difX,
difY;
/* the following are for restarting this routine */
/* dot is just for keeping track of how much has been done */
if(refreshFlag == YES) {
if(dot > 0 && dot < dots)
/* dot >0 is most likely to occur when using
command "cont" */
;
}
else {
dot = 0;
iter = 1;
}
/* dot is just for keeping track of how much has been done */
initTime();
itermax = iterMax;
multiplier = iterMultiplier;/* we need this number globally */
level0 = ++level;
if(level == PROCESS && refreshFlag == NO)
boot_crt();
store(y, y + eqn0, eqn1);
oldX = y[X_coord];
oldY = y[Y_coord];
jump = USTEP; /* this is a fraction of the screen's diagonal
length(currently .001) indicating how big a
jump is allowed so that we do not move more
than 1 pixel */
difX = line_x1 - line_x2;
difY = line_y1 - line_y2;
length = sqrt(difX * difX + difY * difY);
/* this will tally the total length plotted so
far; initially it is the length of the
initial line segment */
diagonal = sqrt((X_upper - X_lower) * (X_upper - X_lower)
+ (Y_upper - Y_lower) * (Y_upper - Y_lower));
boxConstants(); /* this defines constants needed by isInBox()
for checking if y[] is in the box */
oldGoodIt = 0;
goodIterate = 1;
iter = 1;
add = USTEP;
store(fracY, y, dim); /* needed for getAdd() */
add = getAdd(add, fracY);/* one iterate; only place called */
for(; level == level0 && iter <= itermax; iter++) {
goodIterate = iter;
/* first we compute the initial jump size "add" */
if(checkFreq () == YES) {
/* means terminate; checkFreq also handles
occasional transmissions of picture to disk
*/
PRINT "process termination: OUT OF TIME \n");
level = level0 - 1;
continue;
}
scr_rowcol(4, 0);
if(printer >= 2) {
if(iter > 1 && length > 0.&& oldLength > 0.) {
entropy = log(length / oldLength);
PRINT
"length(%d)=length of %dth image of segment=%lf \n",
iter - 1, iter - 1, length * diagonal);
PRINT
"Crude Estimate of topological entropy = log(length(%d)/length(%d))=%5.5lf \n"
,iter - 1, iter - 2, entropy);
}
else
PRINT "length(0)=%lf \n", length);
}
oldLength = length;
length = 0.;
/* MAIN INNER LOOP */
if(refreshFlag == YES)
refreshFlag = NO;
else
frac = 0;/* this occurs always except possibly once */
for(; level == level0 && frac < 1.0;
frac = ((frac + add >= 1.0) ? 1.0 : frac + add)) {
checkForInterrupt();
dot++; /* can be deleted; just for keeping track of
how much has been done and how fast it is */
/* when frac == 0, goodIterate and oldGoodIterate will not be
equal since goodIterate has just been increased beyond
previous values */
getGoodIt();
/* this is the only place in the "for(frac..."
loop that goodIterate is set; this also sets
y according to frac and iterates as many
times up to iter as possible while staying
in the box; at the end, y[] is this final
image point that is in the box and
goodIterate is the number of good iterates
*/
while(goodIterate >= oldGoodIt + 2) {
scr_rowcol(7, 0);
PRINT
"good Iterate just jumped by 2\n");
frac = frac - add;
add = add / 5.;
frac = frac + add;
getGoodIt();
}
store(fracY, y, dim);/* fracY is used by revise... */
if (printer >= 3)
{
scr_rowcol(9,0);
PRINT
"%d good; iter=%d frac=%lf,add=%le \n",goodIterate,iter,frac,add);
}
if (printer >= 1 && minAdd > add)
{
scr_rowcol(8,0);
PRINT
"min add so far =%le \n",add);
}
minAdd = (minAdd > add ? add : minAdd ) ;
reviseAddAndJump(goodIterate, fracY);
/* if(goodIterate <= iter-1) we test to see if frac+add gives
a point that can be iterated one more time; if it can be,
then we will be decrease add by implementing the instruction
: reviseAddAndJump(goodIterate+1,fracY) which should produce
a smaller "add" and just to be safer, we "add = add/3.;" */
if(goodIterate <= iter - 1) {
process(1);
/* one more iterate from y[] as finally set by
revise...() */
if(isInBox (y[X_coord], y[Y_coord]) == YES) {
store(y, fracY, dim);
process(1);
/* one more iterate from y[] as finally set by
GetGoodIt() */
store(fracY, y, dim);
/* now fracY is as set by GetGoodIt() PLUS 1
iterate */
reviseAddAndJump(goodIterate + 1, fracY);
add = add / 3.;
}
}
if(goodIterate == iter || oldGoodIt == iter) {
store(y, fracY, dim);
plotSeg();
length = length + dif;
}
oldGoodIt = goodIterate;
/* the following two coordinates result from iterating y[] iter
times PROVIDED goodIterate = either iter or iter -1; then
they can be used for plotting */
oldX = y[X_coord];
oldY = y[Y_coord];
}
}
level = level0 - 1;
if(level < PROCESS)
MainMenu();
}
double jumpSize() { /* returns the maximum of the horizontal and
vertical fractions of the screen between y
and(oldX,oldY) */
double fabs(), difX, difY;
difX = fabs((y[X_coord] - oldX) * oneOverX);
difY = fabs((y[Y_coord] - oldY) * oneOverY);
if(difX > difY)
return(difX);
return(difY);
}
reviseAddAndJump(its, pY) /* given an estimate "add" of how much frac
should be changed by, try it out using
MakeAJump() which tells us what fraction of
the screen we actually move; if that move is
too big, (that is, if the updated estimate
is less than .8 of add), then do it again
using the new add; the desired jump is
jump*numPixels so we change add by the
factor jump*numPixels/dif. ALERT: if
goodIterate < its at any point, we stop
cycling */
int its;
double *pY;
{
double MakeAJump();
PreviousAdd = add;
for(;;) { /* change "add"; if it decreases, don't
decrease it too much in one step; but keep
re-evaluating while it decreases; decreasing
too much causes a problem when two points
differ primarily because they are being
calculated mod something */
dif = MakeAJump(its, frac + add, pY);
add =.2 * add +.8 * add * numPixels * jump / dif;
/* if add is decreased by more than 20%, go through another time to
refine add */
if(add <= 0.5 * PreviousAdd) {
add = 0.5 * PreviousAdd;
PreviousAdd = add;
}
else
break;
/* checkForInterrupt(); seems possibly unnecessary */
}
if(add > 2.* PreviousAdd) {/* this multiplier should be at least the
reciprocal of multiplier of PreviousAdd */
add = 2.* PreviousAdd;/* don't increase add too fast */
}
}
double MakeAJump(times, fraction, pY)
/* The process is iterated "times" times and
compares how far the new y[] is in the plane
from the vector pY[], and returns that
distance(unless the distance is too small;
then it returns .00...01 ); */
int times;
double fraction,
*pY;
{
double difX,
difY,
ret;
setY(fraction);
process(times);
difX = y[X_coord] - pY[X_coord];
difY = y[Y_coord] - pY[Y_coord];
ret = sqrt(difX * difX + difY * difY) / diagonal;
if(ret <.000001)
ret =.000001;
return(ret);
}
process(times) /* the inner loop (which itself iterates
multipier times) is iterated up to "iter"
times and in fact = "iter" times unless the
trajectory leaves the box, in which case
iterating is terminated. The position y[] at
exit time is its value after "times"
iterates */
int times;
{
int i,
j;
for(i = 1; i <= times; i++) {
for(j = 0; j < multiplier; j++)
(*iteratee) ();
}
}
int isInBox(xval, yval) /* This routine uses the same arithmetic
procedures as plot to see if
(y[X_coord],y[Y_coord]) is within the window
of ScreenSection and if "diameter" is >= 0
and it is not within the window, the routine
checks to see if it is within
"diameters[ScrnSec]" times the width and
within "diameters[ScrnSec] times the height
of the screen. If it is not close enough to
the screen based on this criterion, the
program returns NO, otherwise YES */
double xval,
yval;
{
if((X_upper - xval) * (X_lower - xval) > ddx
|| (Y_upper - yval) * (Y_lower - yval) > ddy) {
return(NO);
}
return(YES);
}
int isInCRT(xval, yval) /* This routine uses the same arithmetic
procedures as plot to see if
(y[X_coord],y[Y_coord]) is within the window
of ScreenSection. If it is not in the
screen, the program returns NO, otherwise
YES */
double xval,
yval;
{
if((X_upper - xval) * (X_lower - xval) > 0
|| (Y_upper - yval) * (Y_lower - yval) > 0) {
return(NO);
}
return(YES);
}
boxConstants() { /* this defines constants needed by isInBox()
for checking if((X_upper-xval)*(X_lower-xval) >
ddx || (Y_upper-yval)*(Y_lower-yval) > ddy )
the constants ddx and ddy must be reset if
the scale or box is changed */
double d,
dd,
xdif,
ydif;
xdif = X_upper - X_lower;
ydif = Y_upper - Y_lower;
oneOverX = 1 / (X_upper - X_lower);
oneOverY = 1 / (Y_upper - Y_lower);
d = diameters[ScrnSec];
dd = d * (1 + d);
ddx = dd * xdif * xdif;
ddy = dd * ydif * ydif;
}
checkForInterrupt() {
#ifndef MAINFRAME
if(ChkKeyStatus () != 0 || cycle > 0) {
Interrupt();
boxConstants();/* this defines constants needed by isInBox()
for checking if y[] is in the box, and
Interrupt can change the required constants
*/
}
#endif
}
plotSeg() { /* plots one or more points depending on
numPixels and modFlag */
double Xp0,
Yp0;
double jumpSize(), js;
int isInCRT(), nPix = numPixels;
int i;
if(goodIterate == iter - 1 || oldGoodIt == iter - 1)
nPix = 10 * nPix;
js = jumpSize();
if(js >.05 + numPixels *.001)
return;
if(cross0flag == 1) { /* this is for preparing to plot a cross */
plotX = y[X_coord];
plotY = y[Y_coord];
}
if(nPix > 1) { /* this segment can screw up if either or both
points were altered by moduloAB() */
Xp0 = y[X_coord];
Yp0 = y[Y_coord];
/* modFlag == YES means that a number has been changed by the
moduloAB() */
for(i = 1; i <= nPix; i++) {
modFlag = NO;
plotX = (Xp0 * i + oldX * (nPix - i)) / nPix;
plotY = (Yp0 * i + oldY * (nPix - i)) / nPix;
if(isInCRT (plotX, plotY) == YES)
plot(plotX, plotY);
}
}
else { /* if numPixels is not greater than 1 */
modFlag = NO;
if(isInCRT (y[X_coord], y[Y_coord]) == YES)
plot(y[X_coord], y[Y_coord]);
}
}
getGoodIt() { /* sets y according to frac and iterates as
many times up to iter as possible while
staying in the box; at the end, y[] is this
final image point and goodIterate is the
number of iterates */
double goodY[EQUATIONS];
int isInBox();
int i,
j;
goodIterate = 0; /* this was = 1 in a buggy earlier version
9/7/88 */
setY(frac);
store(goodY, y, dim);
for(i = 1; i <= iter; i++) {
for(j = 0; j < multiplier; j++)
(*iteratee) ();
if(isInBox (y[X_coord], y[Y_coord]) == YES) {
store(goodY, y, dim);
goodIterate = i;
}
else {
store(y, goodY, dim);
return;
}
}
}
setY(fraction)
double fraction;
{
int i;
for(i = zeroth; i < zeroth + vec_dim; i++)
y[i] = ya[i] + fraction * (yb[i] - ya[i]);
/*** y[X_coord] = line_x1 + fraction*(line_x2 - line_x1);
y[Y_coord] = line_y1 + fraction*(line_y2 - line_y1); ***/
}
double getAdd(add, pY) /* returns estimate for add */
double add,
*pY;
{
int n;
double MakeAJump();
getGoodIt();
/* now estimate add and then refine the estimate */
for(n = 0; n < 2; n++) {
dif = MakeAJump(goodIterate, frac + add, pY);
/* should be about equal to jump*numPixels */
add = add * numPixels * USTEP * diagonal / dif;
}
return(add);
}
#ifdef OBSOLETE
set_line(xshift1, xshift2, yshift1, yshift2)
/* used for drawing lines and unstable
manifolds */
double xshift1,
xshift2,
yshift1,
yshift2;
{
line_x1 = y[eqn0 + X_coord] + xshift1;
line_x2 = y[eqn0 + X_coord] + xshift2;
line_y1 = y[eqn0 + Y_coord] + yshift1;
line_y2 = y[eqn0 + Y_coord] + yshift2;
}
#endif
LineForU(iterMultiplier) /* this routine sets the ends of a line from ya
to yb for computing an unstable manifold
assuming that period NP in Newton's Method
divides the period of the orbit and y1[] is
at a point of the orbit */
int iterMultiplier;
{
double diagonal,
distance(), ratio;
double yy[EQUATIONS],
difX,
difY,
dif;
int i,
J;
diagonal = sqrt((X_upper - X_lower) * (X_upper - X_lower)
+ (Y_upper - Y_lower) * (Y_upper - Y_lower));
/*******
if(printer >= 3)
PRINT
"initial:y1(x) = %11.11lf, y() = %11.11lf, dif=%11.11lf\n"
,y[eqn0+X_coord],y[X_coord],y[X_coord]-y[eqn0+X_coord]);
*********/
/* the following is to set the time coordinates of ya and yb so they will be
equal; other coordinates copied are now irrelevant */
store(ya, y, zeroth);
store(yb, y, zeroth);
/* iterate until y is more than dif/(*USTEP*diagonal) from y1 */
J = 0; /* counts iterates to prevent loop from never
ending */
do {
store(yy, y, dim);
for(i = 0; i < iterMultiplier; i++)
(*iteratee) ();
difX = y[X_coord] - y[eqn0 + X_coord];
difY = y[Y_coord] - y[eqn0 + Y_coord];
dif = sqrt(difX * difX + difY * difY);
ratio = dif / (numPixels * USTEP * diagonal);
/* desired ratio = 1 */
if(printer >= 1 && J == 0 && ratio > 5.)
/* ratio should start out small */
warning1 ();
/********* if(printer >= 3)
{
PRINT
"J = %d: y1(x) = %11.11lf, y() = %11.11lf, dif=%11.11lf\n"
,J,y[eqn0+X_coord],y[X_coord],
y[X_coord]-y[eqn0+X_coord]);
PRINT
"iteration %d: length = %lf length/desired len.= %lf difX=%lf \n"
,J,dif,ratio,difX);
} *********/
} while(++J < 100 && ratio < 1.);
if(J == 100)
PRINT
"Something is wrong. The point may not be unstable \n");
if(printer >= 2)
PRINT "diagonal = %lf ratio = %lf \n", diagonal, ratio);
/* now rescale yy so that is iterates to a point dif/(USTEP*diagonal) from y1*/
if(printer >= 3)
PRINT "before first rescaling, ratio is now %lf\n", ratio);
if(ratio > 0) {
for(i = zeroth; i < lyapzero; i++)
yy[i] = y[eqn0 + i] + (yy[i] - y[eqn0 + i]) / ratio;
store(y, yy, dim);
for(i = 0; i < iterMultiplier; i++)
(*iteratee) ();
difX = y[X_coord] - y[eqn0 + X_coord];
difY = y[Y_coord] - y[eqn0 + Y_coord];
}
dif = sqrt(difX * difX + difY * difY);
ratio = dif / (numPixels * USTEP * diagonal);/* desired ratio = 1 */
if(printer >= 3)
PRINT "after first rescaling, ratio is now %lf\n", ratio);
/* now rescale again in an effort to get closer to the desired length */
for(i = zeroth; i < lyapzero; i++)
yy[i] = y[eqn0 + i] + (yy[i] - y[eqn0 + i]) / ratio;
store(y, yy, dim);
for(i = 0; i < iterMultiplier; i++)
(*iteratee) ();
difX = y[X_coord] - y[eqn0 + X_coord];
difY = y[Y_coord] - y[eqn0 + Y_coord];
dif = sqrt(difX * difX + difY * difY);
dif = sqrt(difX * difX + difY * difY);
ratio = dif / (numPixels * USTEP * diagonal);/* desired ratio = 1 */
if(printer >= 3)
PRINT "final ratio is now %lf; about 1 is target \n", ratio);
/*******
if(printer >= 3)
PRINT
" J = %d:y1(x) = %11.11lf, y(x) = %11.11lf, dif=%11.11lf\n"
,J,y[eqn0+X_coord],y[X_coord],y[X_coord]-y[eqn0+X_coord]);
if(printer >= 2)
PRINT "diagonal = %lf final ratio = %lf \n",diagonal, ratio);
********/
/* now set ends of segment to be ya and yb */
for(i = zeroth; i < zeroth + vec_dim; i++)
ya[i] = yy[i];
for(i = zeroth; i < zeroth + vec_dim; i++)
yb[i] = y[i];
line_x1 = yy[X_coord];
line_x2 = y[X_coord];
line_y1 = yy[Y_coord];
line_y2 = y[Y_coord];
connect2 (y[eqn0 + X_coord], y[eqn0 + Y_coord], line_x1, line_y1);
/* connect2(line_x1,line_y1,line_x2,line_y2);*/
if(printer >= 2) {
PRINT
"Now we compute images of ya-yb = (%lf,%lf)to(%lf,%lf) \n"
,line_x1, line_y1, line_x2, line_y2);
pause(3.);
}
}
LineForB(y1, y2)
double *y1,
*y2;
{
diagonal = sqrt((X_upper - X_lower) * (X_upper - X_lower)
+ (Y_upper - Y_lower) * (Y_upper - Y_lower));
line_x1 = y1[X_coord];
line_x2 = y2[X_coord];
line_y1 = y1[Y_coord];
line_y2 = y2[Y_coord];
stoplines(); /* for stopping drawing connected lines */
#ifdef X11
connectp(line_x1, line_y1);
#else
x_c = line_x1;
y_c = line_y1;
plot(x_c, y_c);
#endif
connectp(line_x2, line_y2);
return;
}
warning1 () { /* for lineForU */
{
PRINT
"/n********************************** WARNING ******************************\n"
);
PRINT
" THE FIRST STEP IS BIG. PERHAPS y1 IS NOT AT A PERIODIC ORBIT OF PERIOD NP\n"
);
PRINT
" PROGRAM WILL ATTEMPT TO CARRY OUT COMMAND\n");
PRINT
"********************************** WARNING ******************************\n\n"
);
}
}